//==============================================================================
// SequenceAlignment.cc
//
//==============================================================================
//.......1.........2.........3.........4.........5.........6.........7.........8
//

#include "SequenceAlignment.hh"

static long limitingCaseCount=0;

// KJD Note: This should be changed to 0,1000000, and 2000000
// after testing.
#define STable 0
#define XTable 10
#define YTable 20

void StartTiming(int &a);
int  EndTiming(int b);

//------------------------------------------------------
// * Set
//------------------------------------------------------
// 
void
SequenceAlignment::SetScoreParams(int match,
						 int mismatch,
						 int gap)
{
	mMatchScore = match;
	mMismatchScore = mismatch;
	mGapScore = gap;
}

//------------------------------------------------------
// * Set
//------------------------------------------------------
// 
void
SequenceAlignment::SetScoreParams(int match,
											 int mismatch,
											 int gap,
											 int extension)
{
	mMatchScore = match;
	mMismatchScore = mismatch;
	mGapScore = gap;
	mGapExtensionScore = extension;
}


//------------------------------------------------------
// * SequenceAlignment
//------------------------------------------------------
// 
SequenceAlignment::SequenceAlignment(Sequence S,Sequence T)
{
  mS.clear();
  mT.clear();
  mS = S;
  mT = T;
  mRows = S.size()+1;
  mCols = T.size()+1;
  SetScoreParams(MATCH_SCORE,MISMATCH_SCORE,GAP_SCORE);
}

//------------------------------------------------------
// * SequenceAlignment
//------------------------------------------------------
// 
SequenceAlignment::SequenceAlignment(string S,string T)
{
	mS.InitFromString(S);
	mT.InitFromString(T);
	mRows = mS.size()+1;
	mCols = mT.size()+1;
	SetScoreParams(MATCH_SCORE,MISMATCH_SCORE,GAP_SCORE);

	mScoreTable.Allocate(mRows,mCols);
	mPtrTable.Allocate(mRows,mCols);
}

//------------------------------------------------------
// * ~SequenceAlignment
//------------------------------------------------------
// 
//SequenceAlignment::~SequenceAlignment(){
//}

//------------------------------------------------------
// * GetGlobalAlignment
//------------------------------------------------------
// Globally aligns two sequences.  
//
// This version has constant gap penalties, so it doesn't distinguish 
// between:
//
//  caggggactacctggaaaataatac
//  ----ggact---tgg------a-a-
//
// and
//
//  caggggactacctggaaaataatac
//  ----ggact---tggaa--------
//
// An affine gap penalty (among others) would prefer the latter. 
//
void 
SequenceAlignment::GetGlobalAlignment(Sequence &AlignS,Sequence &AlignT,
									  int &score)
{
  mScoreTable.Allocate(mRows,mCols);
  mPtrTable.Allocate(mRows,mCols);

  ComputeGlobalAlignmentScore();
  RecoverGlobalAlignment();
  AlignS = mAlignS;
  AlignT = mAlignT;
  score = mScoreTable[mRows-1][mCols-1];
}

//------------------------------------------------------
// * FindWithLocalAlignment
//------------------------------------------------------
// Use local alignment to determine if the target sequence
// is contained in the source sequence.  
//
int
SequenceAlignment::FindWithLocalAlignment(Sequence &Source,Sequence &Target,
														int &score)
{	
  // Ensure that they are initially clear. 
  mS.clear();
  mT.clear();
  mS = Target;
  mT = Source;
  mRows = Target.size()+1;
  mCols = Source.size()+1;
  mAlignStart = -1;
  mMaxScore = 0;

  // Ensure case is same.
  mS.ToLowerCase();
  mT.ToLowerCase();
  
  SetScoreParams(1,-1,-2);
  
  // This is done here since FindWithLocalAlignment is called with the 
  // source and target directly. 
  mScoreTable.Allocate(mRows,mCols);
  mPtrTable.Allocate(mRows,mCols);
  
  //int Start;
  //StartTiming(Start);
  ComputeLocalAlignmentScore();
  //int ElapsedTime = EndTiming(Start);

  //cout<<"Elapsed time="<<ElapsedTime<<endl;

  RecoverLocalAlignment();
	
  score = mMaxScore;
  if (score > Target.size()){
	 cout<<"Score error."<<endl;
  }

  //mScoreTable.DeAllocate();
  //mPtrTable.DeAllocate();

  return(mAlignStart);
}


//------------------------------------------------------
// * GetLocalAlignment
//------------------------------------------------------
// Recovers the best locally aligned subsequence of two 
// sequences.  
//
void 
SequenceAlignment::GetLocalAlignment(Sequence &AlignS,
									 			Sequence &AlignT,int &score)
{
  mScoreTable.Allocate(mRows,mCols);
  mPtrTable.Allocate(mRows,mCols);

  ComputeLocalAlignmentScore();
  RecoverLocalAlignment();
  AlignS = mAlignT;
  AlignT = mAlignS;
}

#if 0
//------------------------------------------------------
// * Score
//------------------------------------------------------ 
//
inline int 
SequenceAlignment::Score(char A,char B)
{
	int score;
	if (A == B){
		score=mMatchScore;
	}else{
		score=mMismatchScore;
	}
	return(score);
}
#endif

//------------------------------------------------------
// * ComputeLocalAlignmentScore
//------------------------------------------------------
// Optimized version of ComputeAlignmentScore.  Clarity is
// sacrificed, a little bit, for the sake of speed.  
//
// KJD: Apparently there is a bug in this somewhere.  
// Occasionally alignments made with this function 
// return with scores greater than the length of the 
// longest sequence.  
//
//
void 
SequenceAlignment::ComputeLocalAlignmentScore4()
{
	register int i,j,im,jm,diagScore,leftScore,topScore;
	register char *pS,*pT;
	register int diagCell,leftCell,topCell,centerCell;
	register int rows,cols,MaxScore=0,MaxRow=0,MaxCol=0;
	const register int GapScore=mGapScore;
	const register int MatchScore = mMatchScore;
	const register int MismatchScore = mMismatchScore;

	// Native character array access is MUCH faster than vector<char>
	// access.  To illustrate, this single efficiency measure dropped
	// the running time from 8 sec to 3 sec, suggesting that 60% of 
	// the inner loop time was spent accessing these characters. 
	pS = new char[mRows];
	pT = new char[mCols];
	for(i = 0;i<mRows;i++){
	  pS[i]=mS[i];
	}

	for(j=0;j<mCols;j++){
	  pT[j]=mT[j];
	}
	
	rows = mRows;
	cols = mCols;

	// Init first two rows to zero for local alignment
	for(i = 0;i<rows;i++)
		mScoreTable[i][0] = 0;
	for(j = 0;j<cols;j++)
		mScoreTable[0][j]=0;

	//------------
	//-------------- Speed critical section ---------------------------
	//
	// Everything in this loop will happen tens of millions of times
	// so it is critical that the total number of operations be minimized.
	// This justifies the expense of copying everything to local register
	// variables, etc.  The utility of the Table class still prompts me 
	// to use it, but that should really be a valarray or something.  
	//
	im = 0;jm=0;  // Previous values of i and j	
	for(i=1;i<rows;i++){
		diagCell = mScoreTable[im][0];
		leftCell = mScoreTable[i][0];
		for(j=1;j<cols;j++){

		  mPtrTable[i][j]=0;

		  if (pS[im]==pT[jm]){
	  	 		diagScore = diagCell+ MatchScore;
	  	 	}else{
	  	 		diagScore = diagCell+ MismatchScore;
	  	 	}
		 
		 	leftScore = leftCell  + GapScore;	 
		 	topCell = mScoreTable[im][j];
		 	topScore =  topCell  + GapScore;
		 
		 	if (0>diagScore){
	 			if (0>topScore){
					if(0>leftScore){
		  				centerCell=0;
					}else{
		  				centerCell=leftScore;
		  				Set(LEFT,mPtrTable[i][j]);
					}
	 			}else{ // a>b a<c 
					if (topScore>leftScore){
		  				centerCell=topScore;
		  				Set(TOP,mPtrTable[i][j]);
					}else{
		  				centerCell=leftScore;
		  				Set(LEFT,mPtrTable[i][j]);
					}
	 			}
  		 	}else{ // a<b
	 			if (diagScore>topScore){
					if (diagScore>leftScore){
		  				centerCell=diagScore;
		  				Set(DIAG,mPtrTable[i][j]);
					}else{
		  				centerCell=leftScore;
		  				Set(LEFT,mPtrTable[i][j]);
					}
	 			}else{ // a<b b<c
					if (topScore>leftScore){
		  				centerCell=topScore;
		  				Set(TOP,mPtrTable[i][j]);
					}else{
		  				centerCell=leftScore;
		  				Set(LEFT,mPtrTable[i][j]);
					}
	 			}
  		 	}
			// Keep track of the location of the maximum score. 
			if (centerCell > MaxScore){
				MaxScore = centerCell;
				MaxRow = i;
				MaxCol = j;
		 	}
		
			mScoreTable[i][j]=centerCell;
			jm=j;
			diagCell = topCell;
			leftCell = centerCell;
		}
		im=i;
	}
	//-------------- End speed critical section -----------------------
	
	mMaxScore = MaxScore;
	mMaxRow = MaxRow;
	mMaxCol = MaxCol;
	delete [] pS;
	delete [] pT;
}


//------------------------------------------------------
// * ComputeLocalAlignmentScore
//------------------------------------------------------
// Optimized version of ComputeAlignmentScore.  Clarity is
// sacrificed, a little bit, for the sake of speed.  
//
//
void 
SequenceAlignment::ComputeLocalAlignmentScore()
{
	register int MaxScore = 0,MaxRow=0,MaxCol=0;
	register int i,j,im,jm,rows,cols,temp,diagScore,leftScore,topScore;
	register char *pS,*pT;
	
	// Native character array access is MUCH faster than vector<char>
	// access.  To illustrate, this single efficiency measure dropped
	// the running time from 8 sec to 3 sec, suggesting that 60% of 
	// the inner loop time was spent accessing these characters. 
	pS = new char[mRows];
	pT = new char[mCols];
	for(i = 0;i<mRows;i++){
	  pS[i]=mS[i];
	}

	for(j=0;j<mCols;j++){
	  pT[j]=mT[j];
	}
	
	rows = mRows;
	cols = mCols;

	// Init first two rows to zero for local alignment
	for(i = 0;i<rows;i++)
		mScoreTable[i][0] = 0;
	for(j = 0;j<cols;j++)
		mScoreTable[0][j]=0;

	//-------------- Speed critical section ---------------------------
	//
	// Everything in this loop will happen tens of millions of times
	// so it is critical that the total number of operations be minimized.
	// This justifies the expense of copying everything to local register
	// variables, etc.  I believe that this can be optimized further by 
	// combining Max and the subsequent Set functions, and by ordering 
	// the comparisons to take advantage of the fact that the score
	// in a local alignment is most everywhere zero.
	//
	for(i=1;i<rows;i++){
	  for(j=1;j<cols;j++){
		  
		 im = i-1;
		 jm = j-1;
		 diagScore = mScoreTable[im][jm]+ Score(pS[im],pT[jm]);
		 leftScore = mScoreTable[i][jm]  + mGapScore;
		 topScore =  mScoreTable[im][j]  + mGapScore;
		 
		 temp = Max(topScore,leftScore,diagScore,0);

		 // KJD BUG HERE??? Note that temp can be zero at times when 
		 // temp != diagScore, leftScore, OR topScore.  What then???
		 // The pointer in the table gets set to nothing in particular. 
		 // What should happen??
	 
		 mScoreTable[i][j]=temp;

		 mPtrTable[i][j]=0;
		 if (temp == topScore) Set(TOP,mPtrTable[i][j]);
		 if (temp == leftScore) Set(LEFT,mPtrTable[i][j]);
		 if (temp == diagScore) Set(DIAG,mPtrTable[i][j]);
		 
		 // Keep track of the location of the maximum score. 
		 if (temp > MaxScore){
			MaxScore = temp;
			MaxRow = i;
			MaxCol = j;
		 }
	  }
	}
	//-------------- End speed critical section -----------------------
	
	mMaxScore = MaxScore;
	mMaxRow = MaxRow;
	mMaxCol = MaxCol;
	delete [] pS;
	delete [] pT;
}


//------------------------------------------------------
// * ComputeLocalAlignmentScore2
//------------------------------------------------------
// This is the original version of ComputeLocalAlignment
// It is retained here for reference as the current 
// ComputeLocalAlignment code evolves in an optimization
// direction, away from comphrehensibility. 
// 
void 
SequenceAlignment::ComputeLocalAlignmentScore2()
{
	register int diagScore,leftScore,topScore,temp;
	register int i,j;

	mScoreTable.Allocate(mRows,mCols);
	mPtrTable.Allocate(mRows,mCols);
	
	// Init first two rows to zero for local alignment
	for(i = 0;i<mRows;i++)
		mScoreTable[i][0] = 0;
	for(j = 0;j<mCols;j++)
		mScoreTable[0][j]=0;

	for(i=1;i<mRows;i++){
	  for(j=1;j<mCols;j++){
			
		 diagScore = mScoreTable[i-1][j-1]+ Score(mS[i-1],mT[j-1]);
		 leftScore = mScoreTable[i][j-1]  + mGapScore;
		 topScore =  mScoreTable[i-1][j]  + mGapScore;

		 temp = Max(diagScore,leftScore,topScore,0);

		 mScoreTable[i][j]=temp;
		 
		 // Save pointers to the places where this value could have
		 // come from. 
		 mPtrTable[i][j]=0;
		 if (temp == topScore) Set(TOP,mPtrTable[i][j]);
		 if (temp == leftScore) Set(LEFT,mPtrTable[i][j]);
		 if (temp == diagScore) Set(DIAG,mPtrTable[i][j]);
		 
		 // Keep track of the location of the maximum score. 
		 if (temp > mMaxScore){
			mMaxScore = temp;
			mMaxRow = i;
			mMaxCol = j;
		 }
	  }
	}
}
		
	
//------------------------------------------------------
// * ComputeGlobalAlignmentScore
//------------------------------------------------------
// 
void 
SequenceAlignment::ComputeGlobalAlignmentScore()
{
	int diagScore,leftScore,topScore,temp;

	mScoreTable.Allocate(mRows,mCols);
	mPtrTable.Allocate(mRows,mCols);
	int i,j;
	for(i = 0;i<mRows;i++)
		mScoreTable[i][0] = i*mGapScore;
	for(j = 0;j<mCols;j++)
		mScoreTable[0][j]=j*mGapScore;

	for(i=1;i<mRows;i++){
		for(j=1;j<mCols;j++){
			
			diagScore = mScoreTable[i-1][j-1]+ 
								Score((char)mS[i-1],(char) mT[j-1]);
			leftScore = mScoreTable[i][j-1]  + mGapScore;
			topScore =  mScoreTable[i-1][j]  + mGapScore;
			
			temp=max(diagScore,leftScore);
			temp=max(temp,topScore);
			mScoreTable[i][j]= temp;
			
			// Save pointers to the places where this value could have
			// come from. 
			if (temp == diagScore) Set(DIAG,mPtrTable[i][j]);
			if (temp == leftScore) Set(LEFT,mPtrTable[i][j]);
			if (temp == topScore) Set(TOP,mPtrTable[i][j]);
		}
	}
}


//------------------------------------------------------
// * Swap
//------------------------------------------------------
void
SequenceAlignment::Swap(int &a,int &b)
{
  int temp = a;
  a = b;
  b = temp;
}


//------------------------------------------------------
// * RecoverGlobalAlignment
//------------------------------------------------------
// Reconstruction of alignment from pointers. 
//  
void 
SequenceAlignment::RecoverGlobalAlignment()
{
  int r = mRows-1,c =mCols-1;
  
  // Reconstruct the reverse, then flip it. 
  while((r >0) && (c >0)){
	 if (FlagSet(DIAG,mPtrTable[r][c])){
		mAlignS.push_back(mS[r-1]);
		mAlignT.push_back(mT[c-1]);
		c--;
		r--;
	 }else{
		if (FlagSet(TOP,mPtrTable[r][c])){
		  mAlignS.push_back(mS[r-1]);
		  mAlignT.push_back(GAP);
		  r--;
		}else{
		  if (FlagSet(LEFT,mPtrTable[r][c])){
			 mAlignS.push_back(GAP);
			 mAlignT.push_back(mT[c-1]);
			 c--;
		  }
		}
	 }
  }
	
  // Since we fill back to front, when strings are 
  // substantially different size the beginning will be
  // filled with blanks.  This fills in the alignment for 
  // this section. (only one of r or c should be > 0). 
  if (r > 0){
	 for(int i = r;i> 0;i--) {
		mAlignS.push_back(mS[i-1]);
		mAlignT.push_back(GAP);
	 }
  }
	
  if (c > 0){
	 for(int i = c;i> 0;i--) {
		mAlignS.push_back(GAP);
		mAlignT.push_back(mT[i-1]);
	 }
  }

  // Reverse, since we constructed a backward version.  
  mAlignS.Reverse();
  mAlignT.Reverse();
}

//------------------------------------------------------
// * RecoverLocalAlignment
//------------------------------------------------------
// Reconstruction of alignment from pointers. 
//  
// KJD: Assumes that larger string is in mT
// 
void 
SequenceAlignment::RecoverLocalAlignment()
{
	int r = mMaxRow,c = mMaxCol;

	mAlignS.clear();
	mAlignT.clear();
 
	// The alignment end is with resptect to the longer string, 
	// here assumed to be the one that corresponds to columns. 
	mAlignEnd = mMaxCol;
 
	while((r > 0) && (c > 0) && (mScoreTable[r][c]!=0)){
	  if (FlagSet(DIAG,mPtrTable[r][c])){
		 mAlignS.push_back(mS[r-1]);
		 mAlignT.push_back(mT[c-1]);
		 c--;
		 r--;
	  }else{
		 if (FlagSet(TOP,mPtrTable[r][c])){
			mAlignS.push_back(mS[r-1]);
			mAlignT.push_back(GAP);
			r--;
		 }else{
			if (FlagSet(LEFT,mPtrTable[r][c])){
			  mAlignS.push_back(GAP);
			  mAlignT.push_back(mT[c-1]);
			  c--;
			}
		 }
	  }		
	}

	//if (mAlignS.size() < mS.size()){
	//  cout<<"Alignment Size Mismatch :"<<r<<","<<c<<endl;
	//  mScoreTable.Print(r-1,r+10,c-1,c+10);
	//}else{
	//  cout<<"Alignment Size OK :"<<r<<","<<c<<endl;
	//  mScoreTable.Print(r-1,r+10,c-1,c+10);
	//}

	// The reverse has been constructed, via push_back, so we 
	// flip it to get the original sense. 
	mAlignS.Reverse();
	mAlignT.Reverse();
	
	// Compute the start of the alignment, which is the end minus
	// the size of the shorter string. 
	mAlignStart = mAlignEnd - mAlignS.size();
}


//------------------------------------------------------
// * PrintSpacePaddedNum
//------------------------------------------------------
//
void
SequenceAlignment::PrintSpacePaddedNum(int x)
{
	// Add a space to all positive numbers to match negative nums. 
	if (x>=0){
		cout<<" ";
	}

	if (abs(x) < 10){
		cout<<"   ";
	}else{
		if (abs(x) < 100){
			cout<<"  ";
		}else{
			if (abs(x) < 1000){
				cout<<" ";
			}
		}
	}
	cout<<x;
}

//------------------------------------------------------
// * PrintScoreTable
//------------------------------------------------------
//
//
void
SequenceAlignment::PrintScoreTable()
{
	int i,j;
	
	// print heading. 
	cout<<"            ";
	for(j = 0;j<mT.size();j++){
		cout<<mT[j]<<"     ";
	}
	cout<<"\n";
	
	for(i = 0;i<mRows;i++){
		if (i > 0){
			cout<<mS[i-1]<<":";
		}else{
			cout<<"  ";
		}
		for(j = 0;j<mCols;j++){
			PrintSpacePaddedNum(mScoreTable[i][j]);
			cout<<" ";
		}
		cout<<"\n";
	}
}

//------------------------------------------------------
// * PrintPtrTable
//------------------------------------------------------
//
void
SequenceAlignment::PrintPtrTable()
{
	int i,j;
	
	// print heading. 
	cout<<"            ";
	for(j = 0;j<mT.size();j++){
		cout<<mT[j]<<"      ";
	}
	cout<<"\n";

	for(i = 0;i<mRows;i++){
		if (i > 0){
			cout<<mS[i-1]<<":";
		}else{
			cout<<"  ";
		}
		cout<<"   ";
		for(j = 0;j<mCols;j++){
			if (FlagSet(LEFT,mPtrTable[i][j])) cout<<"L";
			else cout<<".";
				
			if (FlagSet(DIAG,mPtrTable[i][j])) cout << "D";
			else cout<<".";
			
			if (FlagSet(TOP,mPtrTable[i][j])) cout<<"T   ";
			else cout<<".   ";
			cout<<" ";
		}
		cout<<"\n";
	}
}
	

//------------------------------------------------------
// * Clear
//------------------------------------------------------
// 
void
SequenceAlignment::Clear(int signal,int &c)
{
  c = c & (~signal);
}

//------------------------------------------------------
// * IsSet
//------------------------------------------------------
// 
bool
SequenceAlignment::FlagSet(int signal,int c)
{
  return(c & signal);
}

//------------------------------------------------------
// * IsSet
//------------------------------------------------------
// 
bool
SequenceAlignment::FlagSet(int signal,char c)
{
  return(c & signal);
}


//------------------------------------------------------
// * Set
//------------------------------------------------------
// 
void
SequenceAlignment::Set(int signal,char &c)
{
  c = c | signal;
}

//------------------------------------------------------
// * Set
//------------------------------------------------------
// 
void
SequenceAlignment::Set(int signal,int &c)
{
  c = c | signal;
}

//------------------------------------------------------
// * Print
//------------------------------------------------------
//
// 
void
SequenceAlignment::Print(int start,int end)
{
	int i;
	for(i = start;i<end;i++){
		cout<<mAlignS[i];
	}
	cout<<"\n";
	for(i = start;i<end;i++){
		cout<<mAlignT[i];
	}
	cout<<"\n";
}

//------------------------------------------------------
// * Print
//------------------------------------------------------
// 
void
SequenceAlignment::Print()
{
	int i;

	for(i = 0;i<mAlignT.size();i++){
		cout<<mAlignT[i];
	}
	cout<<"\n";
	for(i = 0;i<mAlignS.size();i++){
		cout<<mAlignS[i];
	}
	cout<<"\n";
}

//------------------------------------------------------
// * Max
//------------------------------------------------------
// 
inline int
SequenceAlignment::Max(int a,int b,int c)
{
  if (a>b){
	 if (a>c){
		return(a);
	 }
  }else{ // a < b
	 if (b>c){
		return(b);
	 }else{ // a<b b<c
		return(c);
	 }
  }
  return(c);
}


//------------------------------------------------------
// * Max
//------------------------------------------------------
// 
inline int
SequenceAlignment::Max(int d,int b,int c,int a)
{
  if (a>b){
	 if (a>c){
		if(a>d){
		  return(a);
		}else{
		  return(d);
		}
	 }else{ // a>b a<c 
		if (c>d){
		  return(c);
		}else{
		  return(d);
		}
	 }
  }else{ // a<b
	 if (b>c){
		if (b>d){
		  return(b);
		}else{
		  return(d);
		}
	 }else{ // a<b b<c
		if (c>d){
		  return(c);
		}else{
		  return(d);
		}
	 }
  }
}

